U okviru ovih vežbi biće predstavljene mogućnosti definisanja i manipulacije nad definicijom referentnog koordinatnog sistema prostornih podataka u R-u. Zatim biće predstavljene mogućnosti učitavanja, obrade, analize i manipulacije nad rasterskim podacima.
Coordinate reference system (CRS) - Koordinatni referentni sistem ili Prostorni referentni sistem - Spatial reference system (SRS) predestavlja sinonim kojim se definiše niz parametara kojim se definiše položaj određene tačke u prostoru u odnosu na položaj Zemlje.
Geografski koordinatni referentni sistem je geometrijski model u kome je definisano sledeće:
Koordinatni referentni sistem u projekciji je geometrijski model u kome je definisano sledeće:
Prema Lott-u 2015 važno je razlikovati:
Učitavanje potrebnih paketa: ***
library(sp)
library(raster)
library(sf)
library(mapview)
library(stringr)
library(tidyverse)
library(dplyr)
library(magrittr)
library(RStoolbox)
library(gridExtra)
library(viridis)
library(stars)## Reading layer `cities' from data source `C:\Users\pbursac\Documents\R\win-library\4.0\rgdal\vectors\cities.shp' using driver `ESRI Shapefile'
## Simple feature collection with 606 features and 4 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -165.27 ymin: -53.15 xmax: 177.1302 ymax: 78.2
## geographic CRS: WGS 84
## Simple feature collection with 606 features and 4 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -165.27 ymin: -53.15 xmax: 177.1302 ymax: 78.2
## geographic CRS: WGS 84
## First 10 features:
## NAME COUNTRY POPULATION CAPITAL geometry
## 1 Murmansk Russia 468000 N POINT (33.08604 68.96355)
## 2 Arkhangelsk Russia 416000 N POINT (40.64616 64.52067)
## 3 Saint Petersburg Russia 5825000 N POINT (30.45333 59.95189)
## 4 Magadan Russia 152000 N POINT (150.78 59.571)
## 5 Perm' Russia 1160000 N POINT (56.23246 58.00024)
## 6 Yekaterinburg Russia 1620000 N POINT (60.61013 56.84654)
## 7 Nizhniy Novgorod Russia 2025000 N POINT (43.94067 56.28968)
## 8 Glasgow UK 1800000 N POINT (-4.269948 55.86281)
## 9 Kazan' Russia 1140000 N POINT (49.14547 55.73301)
## 10 Chelyabinsk Russia 1325000 N POINT (61.39261 55.145)
Kreiaranje data frame-a gde su prostorni podaci - koordinate tačaka date kao dve kolone:
## NAME COUNTRY POPULATION CAPITAL X Y
## 1 Murmansk Russia 468000 N 33.086040 68.96355
## 2 Arkhangelsk Russia 416000 N 40.646160 64.52067
## 3 Saint Petersburg Russia 5825000 N 30.453327 59.95189
## 4 Magadan Russia 152000 N 150.780014 59.57100
## 5 Perm' Russia 1160000 N 56.232464 58.00024
## 6 Yekaterinburg Russia 1620000 N 60.610130 56.84654
## 7 Nizhniy Novgorod Russia 2025000 N 43.940670 56.28968
## 8 Glasgow UK 1800000 N -4.269948 55.86281
## 9 Kazan' Russia 1140000 N 49.145466 55.73301
## 10 Chelyabinsk Russia 1325000 N 61.392612 55.14500
sf paket - definicija CRS-aDefinisanje CRS-a se postiže na sledeći način:
crs_wgs84 <- st_crs(4326) # WGS84 elipsoidni koordinatni sistem ima EPSG kod 4326
class(crs_wgs84) # klasa## [1] "crs"
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
## [1] 4326
## [1] "+proj=longlat +datum=WGS84 +no_defs"
Dodeljivanje CRS-a se postiže na sledeći način:
gradovi2 <- st_as_sf(gradovi_df, coords = c("X", "Y")) # kreiranje objekta sf klase
gradovi2 # CRS: NA## Simple feature collection with 606 features and 4 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -165.27 ymin: -53.15 xmax: 177.1302 ymax: 78.2
## CRS: NA
## First 10 features:
## NAME COUNTRY POPULATION CAPITAL geometry
## 1 Murmansk Russia 468000 N POINT (33.08604 68.96355)
## 2 Arkhangelsk Russia 416000 N POINT (40.64616 64.52067)
## 3 Saint Petersburg Russia 5825000 N POINT (30.45333 59.95189)
## 4 Magadan Russia 152000 N POINT (150.78 59.571)
## 5 Perm' Russia 1160000 N POINT (56.23246 58.00024)
## 6 Yekaterinburg Russia 1620000 N POINT (60.61013 56.84654)
## 7 Nizhniy Novgorod Russia 2025000 N POINT (43.94067 56.28968)
## 8 Glasgow UK 1800000 N POINT (-4.269948 55.86281)
## 9 Kazan' Russia 1140000 N POINT (49.14547 55.73301)
## 10 Chelyabinsk Russia 1325000 N POINT (61.39261 55.145)
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
sp paket - definicija CRS-aDefinisanje CRS-a se postiže na sledeći način:
crs_wgs84 <- CRS(SRS_string = "EPSG:4326") # WGS84 elipsoidni koordinatni sistem ima EPSG kod 4326
class(crs_wgs84) # klasa## [1] "CRS"
## attr(,"package")
## [1] "sp"
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]
Dodeljivanje CRS-a se postiže na sledeći način:
gradovi2 <- gradovi_df # kreiranje objekta sp klase
coordinates(gradovi2) <- ~ X + Y
proj4string(gradovi2) <- crs_wgs84
gradovi2## class : SpatialPointsDataFrame
## features : 606
## extent : -165.27, 177.1302, -53.15, 78.2 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 4
## names : NAME, COUNTRY, POPULATION, CAPITAL
## min values : Abidjan, Afghanistan, -99, N
## max values : Zibo, Zimbabwe, 23620000, Y
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 606 obs. of 4 variables:
## .. ..$ NAME : chr [1:606] "Murmansk" "Arkhangelsk" "Saint Petersburg" "Magadan" ...
## .. ..$ COUNTRY : chr [1:606] "Russia" "Russia" "Russia" "Russia" ...
## .. ..$ POPULATION: num [1:606] 468000 416000 5825000 152000 1160000 ...
## .. ..$ CAPITAL : chr [1:606] "N" "N" "N" "N" ...
## ..@ coords.nrs : int [1:2] 5 6
## ..@ coords : num [1:606, 1:2] 33.1 40.6 30.5 150.8 56.2 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:606] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "X" "Y"
## ..@ bbox : num [1:2, 1:2] -165.3 -53.2 177.1 78.2
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "X" "Y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]
sf <–> sp## Simple feature collection with 606 features and 4 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -165.27 ymin: -53.15 xmax: 177.1302 ymax: 78.2
## geographic CRS: WGS 84 (with axis order normalized for visualization)
## First 10 features:
## NAME COUNTRY POPULATION CAPITAL geometry
## 1 Murmansk Russia 468000 N POINT (33.08604 68.96355)
## 2 Arkhangelsk Russia 416000 N POINT (40.64616 64.52067)
## 3 Saint Petersburg Russia 5825000 N POINT (30.45333 59.95189)
## 4 Magadan Russia 152000 N POINT (150.78 59.571)
## 5 Perm' Russia 1160000 N POINT (56.23246 58.00024)
## 6 Yekaterinburg Russia 1620000 N POINT (60.61013 56.84654)
## 7 Nizhniy Novgorod Russia 2025000 N POINT (43.94067 56.28968)
## 8 Glasgow UK 1800000 N POINT (-4.269948 55.86281)
## 9 Kazan' Russia 1140000 N POINT (49.14547 55.73301)
## 10 Chelyabinsk Russia 1325000 N POINT (61.39261 55.145)
## class : SpatialPointsDataFrame
## features : 606
## extent : -165.27, 177.1302, -53.15, 78.2 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 4
## names : NAME, COUNTRY, POPULATION, CAPITAL
## min values : Abidjan, Afghanistan, -99, N
## max values : Zibo, Zimbabwe, 23620000, Y
raster-aRaster ili grid je prostorna struktura podataka koja deli prostor u pravougaone ćelije (piksele) koje sadrže jednu ili više vrednosti nekog opažanog ili modeliranog fenomena, vrednosti ćelija su uglavnom numerički ili kategorijski podaci.
Više na: Geovizualizacija i Web kartografija, M. Kilibarda i D. Protić - http://osgl.grf.bg.ac.rs/books/gvvk/.
Raster (grid) koji ima jednu atributnu vrednost podseća na matricu koja dodatno ima definisano i zaglavlje. Postavlja se pitanje kako se definiše georeferencija kod rasterskih podataka. Praktično ona je sadržana u zaglavlju. Zavisno od formata rasterskih podataka, u zaglavlju se najčešće nalaze podaci o koordinatnom referentnom sistemu (kod nekih formata su u zasebnom fajlu), broju vrsta i kolonan matrice podataka, koordinatama početne ćelije (piksela), veličini piksela (neki formati podržavaju samo kvadratne dok drugi podržavaju i pravougaone piksele), informacija kako je zapisana vrednost za piksele koji nemaju atributnu vrednost.
U okviru R okruženja, za rad sa prostornim podacima u formi rastera, dostupni su paketi raster i stars.
rasterPaket raster sadrži više klasa za rasterske podatke od kojih su tri glavne RasterLayer, RasterBrick i RasterStack. Raster paket sadrži funkcionalnosti za kreiranje, čitanje, pisanje i manipulaciju nad rasterskim podacima, kao što je na primer rasterska algebra. Isto tako paket sadrži i funkcionalnosti za vizualizaciju rasterskih podataka.
Klase:
RasterLayer predstavlja rasterski podatak sa jednim atributom. Objekat tipa RasterLayer sadrži pored atributnih podataka i sve parametre koje se tipično čuvaju u zaglavlju (hederu) rasterskog fajla. To su: broj vrsta i kolona, prostorni obuhvat podataka i parametri koordinatnog referentnog sistema. Zavisno od veličine fajla svi podaci mogu biti u RAM memoriji ili im se pristupa sa diska računara.RasterBrick je klasa koja podržava više atributa u okviru rasterskih podataka, tipično to važi kad se koristi satelitski snimak sa više kanala, na primer multispektralni snimak iz Sentinel 2 misije ili Landsat 8. Ono što je specifično za RasterBrick je da on može da se referiše samo na jedan fajl koji je sačuvan na disku, na primer GeoTIFF fajl u kome su sačuvani kanali nekog multispektralnog satelitskog snimka.RasterStack je lista rasterskih lejera, vrlo slične strukture kao RasterBrick, ali ne postoji ograničinje da se referiše na jedan GIS fajl, već može da bude kompozicija više fajlova koji su sačuvani na disku ili koji su u RAM memoriji računara, ili kombinacija istih. Praktično se može reći da je to lista objekata tipa RasterLayer. Podrazumeva se da je prostorni obuhvat i rezolucija ista za svaki RasterLayer koji je deo RasterStack objekta.starsPaket stars [https://r-spatial.github.io/stars/] nudi klase i metode za učitavanje, obradu, vizualizaciju i čuvanje podataka u formi vektorskih i rasterskih data cube - kocki podataka. U potpunosti se oslanja na tidyverse familiju paketa i logiku rada sa prostornim podacima koja je implementirana u okviru paketa sf.
Kreiranje rastera sa slucajnim vrednostima celija - piksela
## class : Extent
## xmin : 7440500
## xmax : 7475000
## ymin : 4948000
## ymax : 4976500
Funkcija rnorm - vektor slucajnih brojeva sa normalnom raspodelom
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.037175 -0.675902 -0.000835 -0.000604 0.674081 4.119823
## class : RasterLayer
## dimensions : 285, 345, 98325 (nrow, ncol, ncell)
## resolution : 100, 100 (x, y)
## extent : 7440500, 7475000, 4948000, 4976500 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : layer
## values : -4.037175, 4.119823 (min, max)
Definicija referentnog koordinatnog sistema
podrucje_beograda1 <- podrucje_beograda
podrucje_beograda2 <- podrucje_beograda
podrucje_beograda3 <- podrucje_beograda
podrucje_beograda4 <- podrucje_beogradaDodeljivanje CRS-a se vrsi putem funkcije <- crs(), na vise načina koji daju isti rezultat
crs(podrucje_beograda1) <- 3909
crs(podrucje_beograda2) <- "EPSG:3909"
crs(podrucje_beograda3) <- st_crs(3909)$wkt # a WKT string
crs(podrucje_beograda4) <- CRS(SRS_string = "EPSG:3909") # an sp CRS objectKontrola:
## PROJCRS["MGI 1901 / Balkans zone 7 (with axis order normalized for visualization)",
## BASEGEOGCRS["MGI 1901",
## DATUM["MGI 1901",
## ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",3906]],
## CONVERSION["Balkans zone 7",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",21,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9999,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",7500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",18277]],
## CS[Cartesian,2],
## AXIS["easting (Y)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["northing (X)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]
all.equal(wkt(podrucje_beograda1),
wkt(podrucje_beograda2),
wkt(podrucje_beograda3),
wkt(podrucje_beograda4))## [1] TRUE
Vizualizacija:
## class : RasterLayer
## dimensions : 285, 345, 98325 (nrow, ncol, ncell)
## resolution : 100, 100 (x, y)
## extent : 7440500, 7475000, 4948000, 4976500 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +units=m +no_defs
## source : memory
## names : layer
## values : -4.037175, 4.119823 (min, max)
Dati su Multispektralni satelitski podaci - snimci, satelitske misije Sentinel 2, sa sledećim dostupnim kanalima:
U okviru ovog dela vežbe potrebno je:
## [1] "Data/Multispektralni_snimci/0_2015-08-24_S2.tif"
## [2] "Data/Multispektralni_snimci/0_2015-09-13_S2.tif"
## [3] "Data/Multispektralni_snimci/0_2015-10-23_S2.tif"
Prvi način: svaki pojedinačan snimak kao poseban rasterStack putem for petlje
rasList <- list()
for(i in 1:length(files)){
if(file.size(files[i]) > 0){
rasList[[i]] <- raster::stack(files[i])
} else {
rasList[[i]] <- NULL
}
}
rasList## [[1]]
## class : RasterStack
## dimensions : 28, 57, 1596, 13 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 397599.5, 398169.5, 4987613, 4987893 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs
## names : X0_2015.08.24_S2.1, X0_2015.08.24_S2.2, X0_2015.08.24_S2.3, X0_2015.08.24_S2.4, X0_2015.08.24_S2.5, X0_2015.08.24_S2.6, X0_2015.08.24_S2.7, X0_2015.08.24_S2.8, X0_2015.08.24_S2.9, X0_2015.08.24_S2.10, X0_2015.08.24_S2.11, X0_2015.08.24_S2.12, X0_2015.08.24_S2.13
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535
##
##
## [[2]]
## class : RasterStack
## dimensions : 28, 57, 1596, 13 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 397599.5, 398169.5, 4987613, 4987893 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs
## names : X0_2015.09.13_S2.1, X0_2015.09.13_S2.2, X0_2015.09.13_S2.3, X0_2015.09.13_S2.4, X0_2015.09.13_S2.5, X0_2015.09.13_S2.6, X0_2015.09.13_S2.7, X0_2015.09.13_S2.8, X0_2015.09.13_S2.9, X0_2015.09.13_S2.10, X0_2015.09.13_S2.11, X0_2015.09.13_S2.12, X0_2015.09.13_S2.13
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535
##
##
## [[3]]
## class : RasterStack
## dimensions : 28, 57, 1596, 13 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 397599.5, 398169.5, 4987613, 4987893 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs
## names : X0_2015.10.23_S2.1, X0_2015.10.23_S2.2, X0_2015.10.23_S2.3, X0_2015.10.23_S2.4, X0_2015.10.23_S2.5, X0_2015.10.23_S2.6, X0_2015.10.23_S2.7, X0_2015.10.23_S2.8, X0_2015.10.23_S2.9, X0_2015.10.23_S2.10, X0_2015.10.23_S2.11, X0_2015.10.23_S2.12, X0_2015.10.23_S2.13
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535
Drugi način: svaki pojedinačan snimak kao poseban rasterStack putem lapply funkcije
## [[1]]
## class : RasterStack
## dimensions : 28, 57, 1596, 13 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 397599.5, 398169.5, 4987613, 4987893 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs
## names : X0_2015.08.24_S2.1, X0_2015.08.24_S2.2, X0_2015.08.24_S2.3, X0_2015.08.24_S2.4, X0_2015.08.24_S2.5, X0_2015.08.24_S2.6, X0_2015.08.24_S2.7, X0_2015.08.24_S2.8, X0_2015.08.24_S2.9, X0_2015.08.24_S2.10, X0_2015.08.24_S2.11, X0_2015.08.24_S2.12, X0_2015.08.24_S2.13
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535
##
##
## [[2]]
## class : RasterStack
## dimensions : 28, 57, 1596, 13 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 397599.5, 398169.5, 4987613, 4987893 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs
## names : X0_2015.09.13_S2.1, X0_2015.09.13_S2.2, X0_2015.09.13_S2.3, X0_2015.09.13_S2.4, X0_2015.09.13_S2.5, X0_2015.09.13_S2.6, X0_2015.09.13_S2.7, X0_2015.09.13_S2.8, X0_2015.09.13_S2.9, X0_2015.09.13_S2.10, X0_2015.09.13_S2.11, X0_2015.09.13_S2.12, X0_2015.09.13_S2.13
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535
##
##
## [[3]]
## class : RasterStack
## dimensions : 28, 57, 1596, 13 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 397599.5, 398169.5, 4987613, 4987893 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs
## names : X0_2015.10.23_S2.1, X0_2015.10.23_S2.2, X0_2015.10.23_S2.3, X0_2015.10.23_S2.4, X0_2015.10.23_S2.5, X0_2015.10.23_S2.6, X0_2015.10.23_S2.7, X0_2015.10.23_S2.8, X0_2015.10.23_S2.9, X0_2015.10.23_S2.10, X0_2015.10.23_S2.11, X0_2015.10.23_S2.12, X0_2015.10.23_S2.13
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["UTM zone 34N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",21,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",16034]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]
## PROJCRS["WGS 84 / UTM zone 34N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 34N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",21,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["World - N hemisphere - 18°E to 24°E - by country"],
## BBOX[0,18,84,24]],
## ID["EPSG",32634]]
Funkcija za vizualizaciju podataka u formi kolor kompozita i falš kolor kompozita:
Nazivi snimka - pomoćno:
files.names <- list.files("Data/Multispektralni_snimci/", full.names = F,
pattern = ".tif")
files.names %<>% stringr::str_remove(., pattern = ".tif")
files.names## [1] "0_2015-08-24_S2" "0_2015-09-13_S2" "0_2015-10-23_S2"
plotS2parcel <- function(listS2parc = listS2parc, files.names = files.names, id = 1){
c.plot <- ggRGB(listS2parc[[id]], r = 3, g = 2, b = 1) +
labs(x = "Easting [m]", y = "Northing [m]",
title = "Sentinel 2 satellite image - Color composite",
subtitle = paste("Image: ", files.names[id]),
caption = "UBFCE (2020) \n
Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
theme(line = element_blank(),
panel.background = element_blank())
f.plot <- ggRGB(listS2parc[[id]], r = 7, g = 3, b = 2) +
labs(x = "Easting [m]", y = "Northing [m]",
title = "Sentinel 2 satellite image - Falsh color composite",
subtitle = paste("Image: ", files.names[id]),
caption = "UBFCE (2020) \n
Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
theme(line = element_blank(),
panel.background = element_blank())
gridCF <- grid.arrange(c.plot, f.plot, ncol = 1, nrow = 2)
return(gridCF)
}Vizualizacija:
Pogledati: https://www.indexdatabase.de/
(NIR - Red) / (NIR + Red) Kanali: NIR = B08, i = 7 i RED = BO2, i = 3;
Funkcija za računanje NDVI indeksa nad svim snimcima:
rasListNDVI <- list()
for(i in 1:length(rasStack)){
rasListNDVI[[i]] <- (rasStack[[i]][[7]] - rasStack[[i]][[3]]) / (rasStack[[i]][[7]] + rasStack[[i]][[3]])
}
rasListNDVI## [[1]]
## class : RasterLayer
## dimensions : 28, 57, 1596 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 397599.5, 398169.5, 4987613, 4987893 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 0.1295294, 0.7217312 (min, max)
##
##
## [[2]]
## class : RasterLayer
## dimensions : 28, 57, 1596 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 397599.5, 398169.5, 4987613, 4987893 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 0.1917172, 0.7256516 (min, max)
##
##
## [[3]]
## class : RasterLayer
## dimensions : 28, 57, 1596 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 397599.5, 398169.5, 4987613, 4987893 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 0.1663763, 0.8106579 (min, max)
## layer
## Min. 0.1295294
## 1st Qu. 0.2300365
## Median 0.2941529
## 3rd Qu. 0.3932392
## Max. 0.7217312
## NA's 0.0000000
## layer
## Min. 0.1917172
## 1st Qu. 0.2799979
## Median 0.4163690
## 3rd Qu. 0.5377995
## Max. 0.7256516
## NA's 0.0000000
## layer
## Min. 0.1663763
## 1st Qu. 0.2133507
## Median 0.2371211
## 3rd Qu. 0.3487042
## Max. 0.8106579
## NA's 0.0000000
plotS2parcel <- function(listS2parc = listS2parc, listS2parcNDVI = listS2parcNDVI, files.names = files.names, id = 1){
c.plot <- ggRGB(listS2parc[[id]], r = 3, g = 2, b = 1) +
labs(x = "Easting [m]", y = "Northing [m]",
title = "Sentinel 2 satellite image - Color composite",
subtitle = paste("Image: ", files.names[id]),
caption = "UBFCE (2020) \n
Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
theme(line = element_blank(),
panel.background = element_blank())
f.plot <- ggRGB(listS2parc[[id]], r = 7, g = 3, b = 2) +
labs(x = "Easting [m]", y = "Northing [m]",
title = "Sentinel 2 satellite image - Falsh color composite",
subtitle = paste("Image: ", files.names[id]),
caption = "UBFCE (2020) \n
Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
theme(line = element_blank(),
panel.background = element_blank())
star_pred <- st_as_stars(listS2parcNDVI[[id]])
ndvi.plot <- ggplot()+
geom_stars(data = star_pred)+
scale_fill_viridis(option = "D", na.value = NA, name = "NDVI: ")+
labs(x = "Easting [m]", y = "Northing [m]",
title = "Sentinel 2 satellite image - NDVI",
subtitle = paste("Image: ", files.names[id]),
caption = "UBFCE (2020) \n
Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
theme(line = element_blank(),
panel.background = element_blank(),
legend.position = "bottom")
df <- as.data.frame(listS2parcNDVI[[id]])
hist.plot <- ggplot(data=df, aes(layer)) +
geom_histogram(binwidth = 0.005,
col="red",
aes(y = ..density..,
fill = ..count..)) +
scale_fill_gradient("Count", low = "blue", high = "red")+
stat_function(fun = dnorm,
color = "black",
size = 1.5,
args = list(mean = mean(df$layer), sd = sd(df$layer)))+
labs(x = "NDVI values", y = "Count",
title = "Sentinel 2 satellite image - NDVI histogram",
subtitle = paste("Image: ", files.names[id]),
caption = "UBFCE (2020) \n
Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
theme(line = element_blank(),
panel.background = element_blank(),
legend.position = "bottom")+
annotate("text",x = 0.5, y = 5.5, label = paste("Mean:",round(mean(df$layer),digits = 2)))+
annotate("text",x = 0.5, y = 5, label = paste("SD:",round(sd(df$layer),digits = 2)))
gridCF <- grid.arrange(c.plot, f.plot, ndvi.plot, hist.plot, ncol = 2)
return(gridCF)
}ZADATAK
+Napraviti kratak izveštaj (RMarkodown) sa vizualizacijom i računanjem indeksa hlorofila. +Sačuvati podatke u formatu GeoTIFF za razmenu i čuvanje rasterskih prostornih podataka.